home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / primes.mor < prev    next >
Internet Message Format  |  1995-03-23  |  12KB
  1. From comp.sys.handhelds Thu Jan 31 11:09:44 1991
  2. Path: mentor.cc.purdue.edu!noose.ecn.purdue.edu!samsung!spool.mu.edu!uunet!mcsun!hp4nl!charon!jurjen
  3. From: jurjen@cwi.nl (Jurjen NE Bos)
  4. Newsgroups: comp.sys.handhelds
  5. Subject: HP48: The ultimate factoring program (up to 18 digits)
  6. Summary: A fast program factoring large integers
  7. Keywords: factoring, primes
  8. Message-ID: <2869@charon.cwi.nl>
  9. Date: 31 Jan 91 09:49:50 GMT
  10. Sender: news@cwi.nl
  11. Organization: STORC, Veldhoven
  12. Lines: 250
  13. Originator: jurjen@lijster.cwi.nl
  14.  
  15.  
  16. This is a program that can factor large integers (up to 63 bits) on the
  17. HP48.  It is inspired by the postings of Klaus Kalb of last week.  In fact, I
  18. used two of Klaus' routines (see below).
  19. The program has the following features:
  20. - very fast (average 4 seconds for numbers of 12 digits,
  21.   1 minute for 18 digits.)
  22. - uses two different algorithms for factoring (trial division and Pollard rho)
  23.   and proves all factors prime (using Selfridge's theorem).
  24. - combines RPL (the language in the manual) with forth (what's in the ROM) and
  25.   machine code.  This gives speed and servicability at the same time.
  26. - The recursive calls are nicely shown during the computation.
  27. - A separate prime testing routine.
  28.  
  29. There's also a bug:
  30. - It doesn't quite work with wordsizes that are a multiple of four.  This is
  31.   due to the carries.  If you set the wordsize to 64 bits, it will be
  32.   modified to 63 to make it work.
  33.  
  34. How to use:
  35. Download the object at the end of this posting.  Use ASC-> to convert it into
  36. a directory.  Save it under the name FACTOR.  Enter the directory.  Press
  37. DEMO.  See what happens.
  38. Remark: factoring numbers with large factors (more than 7 digits) can take
  39. up to ten minutes.
  40. Mind the bug above: if you set the wordsize to a multiple of four, factoring
  41. might take forever.
  42.  
  43. Some interesting numbers to try:
  44. 3^37+1 (to get this, type 63 STWS #3 #37 #1 NEG POWMOD 1 +).  The program
  45.   thinks it found a large prime factor, but quickly finds out it isn't.
  46. 100264053529.  This number looks even more prime.
  47.  
  48. (Intermezzo:)
  49. Some puzzles for experienced HP48 users:
  50.     (For information, send Email.  Please realize these are PUZZLES, so
  51.     they aren't trivial, and explanation is complicated.)
  52. - Compute 'SIN(180_\^o)' .  You do not get 0.  Why?
  53. - The time to compute '253!' is much longer than the time to compute '253.1!'.
  54. - To compute the factorial for large negative x, use '\pi*x/(SIN(\pi*x)*(-x)!)'
  55.   instead of 'x!' for faster results.  (Don't forget RAD.)
  56.  
  57. Happy Hacking,
  58.     Jurjen Bos
  59.  
  60. -------------------Listing of the FACTOR directory (do not download)------
  61. For the nosy people among you, here the listing of the routines:
  62. (Don't bother downloading these, it's all in the ASC string at the end.)
  63.  
  64. DEMO    @ Factor a random integer and time the computation
  65. \<< DEC TICKS RAND 1000000 * R\->B RAND 1,E12 * + RAND 1,E18 * + FACTOR
  66.   TICKS ROT - B\->R 1,220703125E-4_s * SWAP
  67. \>>
  68.  
  69. PRIME?    @Prime test
  70. \<< # 0d +
  71.   CASE DUP # 2d <
  72.     THEN DROP 0
  73.     END DUP # 210d BGCD # 1d ==
  74.     THEN RCWS 63 MIN STWS 0 'L' STO LPRT 'L' PURGE
  75.     END DUP # 121d <
  76.     THEN B\->R { 2 3 5 7 } SWAP POS 0 \=/
  77.     END DROP 0
  78.   END
  79. \>>
  80.  
  81. FACTOR    @The factoring routine
  82. \<< RCWS 63 MIN STWS # 0d + 0 'L' STO LFAC 'L' PURGE
  83. \>>
  84.  
  85. MULMOD    @Modulo multiplication of binaries
  86. (In machine code.
  87. a b n -> a*b mod n
  88. )
  89.  
  90. POWMOD    @Modulo powering of binaries
  91. (In machine code.
  92. a e n -> a^e mod n
  93. )
  94.  
  95. BGCD    @Greatest commond divisor of binaries
  96. (In machine code.  Written by Klaus Kalb.
  97. )
  98.  
  99. TRIAL    @Trial division routine
  100. (In machine code and forth.  Machine code part written by Klaus Kalb.
  101. )
  102.  
  103. MILRAB    @Miller-Rabin test
  104. (In forth.
  105. n base -> 1 if n is pseudoprime 
  106. )
  107.  
  108. LFAC    @Internal factoring routine
  109. \<< DUP 'L' INCR DISP "Trial division" 'L' INCR DISP # 2000d TRIAL
  110.   IF DUP # 1d ==
  111.   THEN DROP { }
  112.   ELSE 1 \->LIST    @List of composite factors
  113.   END
  114.   WHILE DUP SIZE
  115.   REPEAT SWAP OVER 1 GET
  116.     CASE DUP # 4012009d <
  117.       THEN B\->R +
  118.       END DUP L 1 - DISP "Primetest" L DISP DUP LPRT
  119.       THEN
  120.         IF DUP # 1000000000000d <
  121.         THEN B\->R
  122.         END +
  123.       END "Pollard \Gr" L DISP
  124.       WHILE P\Gr DUP # 1d ==
  125.       REPEAT DROP
  126.       END 2 \->LIST ROT SWAP + SWAP
  127.     END SWAP 2 OVER SIZE SUB
  128.   END DROP 'L' "
  129. " OVER DECR DISP 1 STO-
  130. \>>
  131.  
  132. LPRT    @Internal primtesting routine
  133. \<<
  134.   CASE DUP # 3d MILRAB NOT
  135.     THEN 0
  136.     END DUP # 25000000000d >
  137.     THEN SRT DUP
  138.     END DUP # 2d MILRAB NOT
  139.     THEN 0
  140.     END DUP # 1373653d <
  141.     THEN 1
  142.     END DUP # 5d MILRAB NOT
  143.     THEN 0
  144.     END DUP # 25326001d <
  145.     THEN 1
  146.     END DUP # 7d MILRAB NOT
  147.     THEN 0
  148.     END DUP # 3215031751d \=/
  149.   END SWAP DROP
  150. \>>
  151.  
  152. SRT    @Test primality of large numbers with Selfridge test.
  153. \<< DUP 1 - LFAC # 3d \-> n l a
  154.   \<< "Selfridge" L DISP 9 CF 1 1 l SIZE
  155.     FOR k
  156.       IF l k GET SWAP OVER \=/ 8 FC?C OR
  157.       THEN
  158.         CASE a n 3 PICK / n POWMOD # 1d DUP2 \=/
  159.           THEN SWAP 3 PICK # 0d + n POWMOD \=/ 8 + SF
  160.           END + 'a' STO+ RAND ,1 \>=
  161.           THEN
  162.           END n a MILRAB
  163.           THEN
  164.           END 9 SF
  165.         END
  166.       END 8 FS? 9 FS? -
  167.     STEP DROP 9 FC?C
  168.   \>>
  169. \>>
  170.  
  171. P\Gr    @Pollard rho factoring algorithm.
  172. \<< \-> n
  173.   \<< RAND n B\->R * R\->B DUP NEWOB
  174.     WHILE n Code n BGCD DUP # 1d ==
  175.     @ Code is a machine code routine that does 50 rounds
  176.     REPEAT DROP
  177.     END ROT ROT DROP2 n
  178.   \>> OVER /
  179. \>>
  180.  
  181. -------------------FACTOR directory in downloadable form------------------
  182. This is the entire program. Download in your machine, use ASC->,
  183. and save under the name FACTOR.  You'll need about 7K free to download it,
  184. it is about 2430 bytes.
  185.  
  186. %%HP: T(3)A(D)F(,);
  187. "69A20FF77E1100000020057920D9D20E16321C432D6E2010E6E16329B1C1D6E2
  188. 010E6BB691EEDA1B969178BF1CB2A133032D6E2010E6CCD208E0008FD8F357C8
  189. 0AF2E610A156110815711091923A98118AF77380B16A977970B16108119A9779
  190. 60B1610911099250A9EB12A9711A7B4010A18013A13896CFA781011015011111
  191. 5111128DD69508F2D76014713416917414713517901A9082281F832D0A1A9945
  192. 0B10A1699150B1991FCDA9601D6E2010E684E20402474344478BF1E4A2060000
  193. 1279E1D50328DBF149632E0CF1E0CF13FBF1D6E2010E6EF53292CF150FA19363
  194. 2B21304B1003035254530D9D20E163278BF19C2A290DA184E2040C4641434E4A
  195. 206000031C432D6E2010E6D6E2010C6D6E201016E1632C2A20710003556C6662
  196. 79646765684E2010C4485A1173A25D2C19C2A29C2A2D6E2010C68B9C10A132D6
  197. E2010B63CE22D6E2010C6D6E2010B66C7D1DBBF192CF1D9AE1C53A2025C1908E
  198. 1AFE22D9D20D8732D9D20D6E201016D6E2010E63F2A2A9CF150FA1D6E2010E68
  199. 4E206005F475D4F444E4A205100010000000000000002ABF1D9AE18A732D9D20
  200. DBBF13F2A2A9CF1E4A2051000000000000000000076BA1D6E2010E684E206005
  201. F475D4F444D9AE1C53A276BA1472C1B21305DF2276BA145632D6E20101697632
  202. B44029B1C1339209990000000000010B9DE18A732D9D20B21305DF22D6E2010E
  203. 6D6E20101684E2060D494C42514248A732D9D20B21305DF22173A2472C1B2130
  204. 5DF22B21305DF22C53A2313C1173A2313C190DA1083328DBF1173A2025C1EF53
  205. 293632B21300C20040C405254540D9D20E1632D8732D9D2078BF1E4A20510003
  206. 00000000000000084E2060D494C4251424F88E18A7324B2A25DF2278BF1E4A20
  207. 5100000ABD12D50000000D5CE18A732D9D2084E203035254578BF1B21305DF22
  208. 78BF1E4A2051000200000000000000084E2060D494C4251424F88E18A7324B2A
  209. 25DF2278BF1E4A20510005D5F410000000000EBBE18A7329C2A25DF2278BF1E4
  210. A2051000500000000000000084E2060D494C4251424F88E18A7324B2A25DF227
  211. 8BF1E4A20510001B17281000000000EBBE18A7329C2A25DF2278BF1E4A205100
  212. 0700000000000000084E2060D494C4251424F88E18A7324B2A25DF2278BF1E4A
  213. 20510007CD71AFB00000000D9AE1B21305DF22DBBF18DBF193632B2130A22004
  214. 0C464143440D9D20E163278BF14563284E2010C4976324F802485A1C2A201200
  215. 045279616C602469667963796F6E64563284E2010C4976324F802485A1E4A205
  216. 10000D7000000000000084E205045259414C43CE2278BF1E4A20510001000000
  217. 000000000279E1AFE22D9D208DBF147A20B2130B21305BF22D9D209C2A2387C1
  218. B21305DF223303278BF18B9C1D5032D9D20DBBF192CF19C2A26C7D1D8732D9D2
  219. 078BF1E4A20510009E73D30000000000EBBE18A732D9D20BB69176BA1B21305D
  220. F2278BF184E2010C49C2A290DA1485A1C2A2071000052796D6564756374784E2
  221. 010C4485A178BF184E2040C40525458A732D9D203CE2278BF1E4A20510000001
  222. 5A4D8E000000EBBE1AFE22BB6915DF2276BA1B21305DF22C2A207100005F6C6C
  223. 6162746027984E2010C4485A13303284E2020057978BF1E4A205100010000000
  224. 00000000279E1D50328DBF149632ED2A2387C1E0CF1DBBF176BA1DBBF1B21305
  225. DF22DBBF1ED2A292CF18B9C1C58C1B2130496328DBF14563284E2010C497632C
  226. 2A2070000A092CF1AA902485A19C2A28350293632B2130F230060D494C425142
  227. 460D9D20FDE8111920BB000D9D202C230E4A206000010BE3588130CCD2042000
  228. 8FD8F3582281C832AFA14B148DD69505AC26A321684E206005F475D4F4448813
  229. 0E4A206000019D445FC7A239916D9D20E7F069C2A2B21302A17088130C12169D
  230. 445B67A2EF116A32169D445B67A264B30EE170D9D2088130A321684E2060D455
  231. C4D4F44432230E5D3532230B21305E17044230CE445B9F06B2130B2130741005
  232. 045259414C450D9D20FDE8111920BB000D9D20CCD20B61008F77F3510110AAF2
  233. 10810B2081B58082444000C21366570242462642466264264684242486462462
  234. 664246264242A2A021224011BD2BF6BF6BF6BF6BF6BF61088F2D7608F735608F
  235. B97601112F8DD6950AF015A097CB12081B580824C7FFFC21366EDF118C24A910
  236. 81181129F2C8111118AF3AF19F262A76B779FE7F81EA75A7F9F280B7AB7597F9
  237. E11A9F180AF910A97CD3AF4101113132AF226301A7A103208F2D7608F735608F
  238. B9760113130657F1606E3F47A20D6E2010B6D6E201027B21300D470D6E2010B6
  239. 6AC30A2170D9D20D6E2010B63C370FBD81D6E2010B65233043370B2130D6E201
  240. 0B695450D6E20102779470B2130B213012200402474344440D9D20FDE81ECD46
  241. CCD20F70008F77F35AF397845AFE978C48087091AFE80870F081E81CB7765EF8
  242. 08B09081E65FF9F650AFEB7A9780181C808608F65EFAFA97BC0A74A7F64FF8F5
  243. 9E35B21304A0006005F475D4F44460D9D200FE8111920BBB00CCD20D80008FD8
  244. F35101208F2D7608F77F35121A98A90B6482281E108832901197C10121A96721
  245. 012111891E8D8DD6950A97A96A9082281F832D0A1A99450B10A1699150B1991F
  246. CD01B2130BB00060D455C4D4F44460D9D200FE8111920BBB00CCD20E50008FD8
  247. F35101208F2D7608F77F35A97119A9182281F83201A1847099550B11A1447099
  248. 250B1A91F6DA948DD6950B2130C80006064143445F42560D9D20E1632EF5C133
  249. 92010000000000003603ECB15C5C1E4A2051000000000000000000076BA14B2A
  250. 24563284E2010C497632DCC0284E2040C46414344563284E2010C497632EFE02
  251. 93632B2130BA00060052594D454F360D9D20E1632E4A20510000000000000000
  252. 00076BA1D8732D9D2078BF1E4A20510002000000000000000EBBE18A732D9D20
  253. 8DBF14B2A2B21305DF2278BF1E4A20510002D0000000000000084E2040247434
  254. 44E4A20510001000000000000000279E18A732D9D20EF5C13392010000000000
  255. 003603ECB15C5C14B2A24563284E2010C497632DCC0284E2040C405254545632
  256. 84E2010C497632EFE02B21305DF2278BF1E4A20510009700000000000000EBBE
  257. 18A732D9D20BB69147A20ED2A23F2A2D13A2743A2B2130DBBF14BAC14B2A2D9A
  258. E1B21305DF228DBF14B2A2B21305DF2293632B2130BD100404454D4F440D9D20
  259. E1632475C1D28919B1C1339206000000000000010EEDA1B96919B1C133920210
  260. 0000000000010EEDA176BA19B1C1339208100000000000010EEDA176BA184E20
  261. 6064143445F425D2891E0CF190DA1BB691ADA20339206990052130702210C2A2
  262. 0700003768B01B2130EEDA1DBBF193632B21306F87"
  263.  
  264. --- end of posting ---
  265.  
  266. From comp.sys.handhelds Thu Jan 31 11:10:04 1991
  267. Path: mentor.cc.purdue.edu!noose.ecn.purdue.edu!samsung!uunet!mcsun!hp4nl!charon!jurjen
  268. From: jurjen@cwi.nl (Jurjen NE Bos)
  269. Newsgroups: comp.sys.handhelds
  270. Subject: HP48: factoring program, patch
  271. Summary: Minor patch.
  272. Keywords: factoring, prime, patch
  273. Message-ID: <2870@charon.cwi.nl>
  274. Date: 31 Jan 91 10:44:12 GMT
  275. Sender: news@cwi.nl
  276. Organization: STORC, Veldhoven.
  277. Lines: 23
  278. Originator: jurjen@lijster.cwi.nl
  279.  
  280.  
  281. The bug mentioned in the previous posting can be succesfully solved by the
  282. following minor patches.  This is typed by hand.  The easiest way to patch is
  283. just using VISIT instead of trying to download it.
  284.  
  285. Happy hacking,
  286.     Jurjen Bos
  287.  
  288. PRIME?
  289. \<< # 0d + # FFFFFFFFFFFFFFFFh SR AND
  290.   CASE DUP # 2d <
  291.     THEN DROP 0
  292.     END DUP # 210d BGCD # 1d ==
  293.     THEN 0 'L' STO LPRT 'L' PURGE
  294.     END DUP # 121d <
  295.     THEN B\->R { 2 3 5 7 } SWAP POS  0 \=/
  296.     END DROP 0
  297.   END
  298. \>>
  299.  
  300. FACTOR
  301. \<< # 0d + # FFFFFFFFFFFFFFFFh SR AND 0 'L' STO LFAC 'L' PURGE
  302. \>>
  303.  
  304.